{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1.4 Static Condensation\n",
    "\n",
    "Static condensation refers to the process of eliminating unknowns that are internal to elements from the global linear system. They are useful in standard methods and critical in methods like the HDG method. NGSolve automates static condensation across a variety of methods via a classification of degrees of freedom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import netgen.gui\n",
    "%gui tk\n",
    "from ngsolve import *\n",
    "from netgen.geom2d import unit_square\n",
    "\n",
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.4))\n",
    "\n",
    "fes = H1(mesh, order=4, dirichlet='bottom|right')\n",
    "u, v = fes.TnT()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Asking BilinearForm to condense"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = BilinearForm(fes, eliminate_internal=True)\n",
    "a += SymbolicBFI (grad(u) * grad(v))\n",
    "a.Assemble()\n",
    "\n",
    "f = LinearForm(fes)\n",
    "f += SymbolicLFI (1 * v)\n",
    "f.Assemble()\n",
    "\n",
    "u = GridFunction(fes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The assembled matrix $A=$ `a.mat` can be block partitioned into\n",
    "\n",
    "$$\n",
    "A =\n",
    "\\left(\n",
    "\\begin{array}{cc}\n",
    "A_{LL} & A_{LE}\n",
    "\\\\\n",
    "A_{EL} & A_{EE}\n",
    "\\end{array}\n",
    "\\right)\n",
    "$$\n",
    "\n",
    "where\n",
    "$L$ denotes the set of **local or internal**  degrees of freedom and $E$ denotes the set of **interface** degrees of freedom.\n",
    "\n",
    "* In our current example $E$ consists of edge and vertex dofs while $L$ consists of triangle dofs. (Note that in practice,  $L$ and $E$ may not be ordered contiguously and $L$ need not appear before $E$, but such details are immaterial for our discussion here.)\n",
    "\n",
    "* The condensed system is known as the **Schur complement**:\n",
    "\n",
    "$$\n",
    "S = A_{EL} - A_{EL} A_{LL}^{-1} A_{LE}.\n",
    "$$\n",
    "\n",
    "* When `eliminate_internal` is set to `True` in `a`, the statement `a.Assemble` actually assembles $S$.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A factorization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "* NGSolve provides \n",
    "\n",
    "    - `a.harmonic_extension_trans` $=\n",
    "\\left(\\begin{array}{cc}\n",
    "0 & 0 \\\\\n",
    "-A_{EI} A_{LL}^{-1} & 0 \n",
    "\\end{array}\\right)$\n",
    "    - `a.harmonic_extension` $=\n",
    "\\left(\\begin{array}{cc}\n",
    "0 & -A_{LL}^{-1} A_{LE}\n",
    "\\\\\n",
    "0 & 0\n",
    "\\end{array}\\right)$\n",
    "\n",
    "    - `a.inner_solve`  $=\\left(\\begin{array}{cc}\n",
    "A_{LL}^{-1} & 0 \n",
    "\\\\\n",
    "0 & 0\n",
    "\\end{array}\\right)$.\n",
    "\n",
    "* To solve $$ \\left(\\begin{array}{cc}\n",
    "A_{LL} & A_{LE}\\\\\n",
    "A_{EL} & A_{EE}\n",
    "\\end{array}\\right)\n",
    "\\left(\\begin{array}{c}\n",
    "u_L \\\\ u_E\n",
    "\\end{array}\\right)= \n",
    "\\left(\\begin{array}{c}\n",
    "f_L \\\\ f_E\n",
    "\\end{array}\\right)\n",
    "$$\n",
    "\n",
    "we use a factorization of $A^{-1}$ that uses $S^{-1}$. Namely, we use the following identity:\n",
    "\n",
    "$$\n",
    "\\left(\\begin{array}{c}\n",
    "u_L \\\\ u_E\n",
    "\\end{array}\\right)\n",
    "=\n",
    "\\left(\\begin{array}{cc}\n",
    "I & -A_{LL}^{-1} A_{LE} \\\\\n",
    "0 & I \n",
    "\\end{array}\\right)\n",
    "\\left(\\begin{array}{cc}\n",
    "A_{LL}^{-1} & 0 \\\\\n",
    "0 & S^{-1}\n",
    "\\end{array}\\right)\n",
    "\\underbrace{\n",
    "\\left(\\begin{array}{cc}\n",
    "I & 0 \\\\\n",
    "-A_{EL} A_{LL}^{-1} & I \n",
    "\\end{array}\\right)\n",
    "\\left(\\begin{array}{c}\n",
    "f_L \\\\ f_E\n",
    "\\end{array}\\right)}_{\n",
    "\\left(\\begin{array}{c}\n",
    "f'_L \\\\ f'_E\n",
    "\\end{array}\\right)\n",
    "}\n",
    "$$\n",
    "\n",
    "We implement this formula step by step, starting with the computation of $f_L'$ and $f_E'$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Steps to compute the solution\n",
    "* The following step implements\n",
    "$$\n",
    "\\left(\\begin{array}{c}\n",
    "f'_L \\\\ f'_E\n",
    "\\end{array}\\right) =\n",
    "\\left(\\begin{array}{cc}\n",
    "I & 0 \\\\\n",
    "-A_{EL} A_{LL}^{-1} & I \n",
    "\\end{array}\\right)\n",
    "\\left(\\begin{array}{c}\n",
    "f_L \\\\ f_E\n",
    "\\end{array}\\right).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "f.vec.data += a.harmonic_extension_trans * f.vec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The next step implements part of the next matrix application in the formula.\n",
    "$$\n",
    "\\left(\\begin{array}{c}\n",
    "0 \\\\ u_E\n",
    "\\end{array}\\right) = \n",
    "\\left(\\begin{array}{cc}\n",
    "0 & 0 \\\\\n",
    "0 & S^{-1}\n",
    "\\end{array}\\right)\n",
    "\\left(\\begin{array}{c}\n",
    "f'_L \\\\ f'_E\n",
    "\\end{array}\\right).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.vec.data = \\\n",
    "a.mat.Inverse(freedofs=fes.FreeDofs(coupling=True)) * f.vec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Note: \n",
    "    - Because we set `eliminate_internal` in `a`, the inverse `a.mat.Inverse` actually computes $S^{-1}$.\n",
    "\n",
    "    - Note that instead of the usual `fes.FreeDofs()`, we have used  `fes.FreeDofs(coupling=True)` or simply `fes.FreeDofs(True)` to specify that only the degrees of freedom that are *not local* and *not Dirichlet* should participate in the inverse computations.  (The underlying assumption is that Dirichlet dofs cannot be local dofs.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Next, we compute\n",
    "$$\n",
    "\\left(\\begin{array}{c}\n",
    "u'_L \\\\ u_E\n",
    "\\end{array}\\right) = \n",
    "\\left(\\begin{array}{c}\n",
    "0 \\\\ u_E\n",
    "\\end{array}\\right) + \n",
    "\\left(\\begin{array}{cc}\n",
    "A_{LL}^{-1} & 0 \\\\\n",
    "0 & 0\n",
    "\\end{array}\\right)\n",
    "\\left(\\begin{array}{c}\n",
    "f'_L \\\\ f'_E\n",
    "\\end{array}\\right).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.vec.data += a.inner_solve * f.vec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Finally:\n",
    "$$\n",
    "\\left(\\begin{array}{c}\n",
    "u_L \\\\ u_E\n",
    "\\end{array}\\right) =\n",
    "\\left(\\begin{array}{cc}\n",
    "I & -A_{LL}^{-1} A_{LE} \\\\\n",
    "0 & I \n",
    "\\end{array}\\right)\n",
    "\\left(\\begin{array}{c}\n",
    "u_L' \\\\ u_E\n",
    "\\end{array}\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.vec.data += a.harmonic_extension * u.vec\n",
    "Draw(u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Behind the scenes: `CouplingType`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does NGSolve know what is in the index sets $L$ and $E$? \n",
    "\n",
    "* Look at `fes.CouplingType` to see a classification of degrees of freedom.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "COUPLING_TYPE.WIREBASKET_DOF\n",
      "COUPLING_TYPE.WIREBASKET_DOF\n",
      "COUPLING_TYPE.WIREBASKET_DOF\n",
      "COUPLING_TYPE.WIREBASKET_DOF\n",
      "COUPLING_TYPE.WIREBASKET_DOF\n",
      "COUPLING_TYPE.WIREBASKET_DOF\n",
      "COUPLING_TYPE.WIREBASKET_DOF\n",
      "COUPLING_TYPE.WIREBASKET_DOF\n",
      "COUPLING_TYPE.WIREBASKET_DOF\n",
      "COUPLING_TYPE.WIREBASKET_DOF\n",
      "COUPLING_TYPE.WIREBASKET_DOF\n",
      "COUPLING_TYPE.WIREBASKET_DOF\n",
      "COUPLING_TYPE.WIREBASKET_DOF\n",
      "COUPLING_TYPE.WIREBASKET_DOF\n",
      "COUPLING_TYPE.WIREBASKET_DOF\n",
      "COUPLING_TYPE.WIREBASKET_DOF\n",
      "COUPLING_TYPE.WIREBASKET_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.INTERFACE_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n",
      "COUPLING_TYPE.LOCAL_DOF\n"
     ]
    }
   ],
   "source": [
    "for i in range(fes.ndof):\n",
    "    print(fes.CouplingType(i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{COUPLING_TYPE.WIREBASKET_DOF: 17,\n",
       " COUPLING_TYPE.INTERFACE_DOF: 108,\n",
       " COUPLING_TYPE.LOCAL_DOF: 60}"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dof_types = {}\n",
    "for i in range(fes.ndof):\n",
    "    ctype = fes.CouplingType(i)\n",
    "    if ctype in dof_types.keys():\n",
    "        dof_types[ctype] += 1\n",
    "    else:\n",
    "        dof_types[ctype] = 1\n",
    "dof_types"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `LOCAL_DOF` forms the set $L$ and the remainder forms the set $E$. All finite element spaces in NGSolve have such dof classification.\n",
    "\n",
    "Through this classification a bilinear form is able to automatically compute the Schur complement and the accompanying extension operators. Users need only specify the  `eliminate_internal` option. (Of course users should also make sure their method has an invertible $A_{LL}$!)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inhomogeneous Dirichlet boundary conditions\n",
    "\n",
    "In case of inhomogeneous Dirichlet boundary conditions we must combine the technique of Dirichlet data extension in a previous tutorial with static condensation:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0192514570333646e-15"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U = x*x*(1-y)*(1-y)          # U = manufactured solution\n",
    "DeltaU = 2*((1-y)*(1-y)+x*x) # Source: DeltaU = ∆U\n",
    "\n",
    "f = LinearForm(fes)\n",
    "f += SymbolicLFI(-DeltaU*v)\n",
    "f.Assemble()\n",
    "\n",
    "u = GridFunction(fes)\n",
    "u.Set(U, BND)               # Dirichlet b.c: u = U\n",
    "\n",
    "# Modify source per Dirichlet extension technique:\n",
    "r = f.vec.CreateVector()\n",
    "r.data = f.vec - a.mat * u.vec\n",
    "r.data += a.harmonic_extension_trans * r\n",
    "\n",
    "# Apply the static condensation technique:\n",
    "u.vec.data += a.mat.Inverse(fes.FreeDofs(True)) * r\n",
    "u.vec.data += a.harmonic_extension * u.vec\n",
    "u.vec.data += a.inner_solve * r\n",
    "\n",
    "sqrt(Integrate((U-u)*(U-u),mesh))  # Compute error"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
